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v .  SUMMARY  „  J  , 
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This  Paper  traces  the  development,  at  the  -ftfe,  of  computer  programs  for 
orbit  determination  and  orbit  evolution.  The  principal  program,  PROP,  is  still 
in  regular  use,  having  an  accuracy  of  50  to  100  m  for  the  orbits  of  close 
satellites  determined  over  a  period  of  a  few  days  even  when  the  drag  is  severe. 
The  chief  limitation  of  PROP  is  the  absence  of  lunisolar  perturbations,  but  these 
are  included  in  a  companion  program,  PROD,  for  the  evolution  of  mean  elements. 

Work  of  recent  years  has  focused  on  a  framework  for  far  more  accurate 
representation  and  generation  of  orbits.  The  method  is  purely  analytical  (as  in 
PROP)  for  short-period  perturbations,  applied  compactly  via  the  use  of  a  special 
system  of  spherical  coordinates,  but  employs  a  degree  of  numerical  integration 
in  the  variation  of  mean  elements.  The  goal  is  a  single  program  that  would 
eventually  supersede  the  PROP-PROD  combination,  providing  an  accuracy  of  perhaps 
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This  Paper  was  presented  (5  July  1984)  at  the  25th  COSPM  Meeting,  held  in 
Graz  (Austria).  It  will  be  published  in  the  COSPAR  Proceedings  (Ado.  Space  Res., 
Vol  4).  The  only  change  from  the  presented  paper  is  the  footnote  on  page  3. 
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EARLY  WORK 


The  Royal  Aircraft  Establishment  became  involved  in  satellite  orbit  determination 
within  a  few  days  of  the  launch  of  Sputnik  1  on  4  October  1957.  Kinetheodolites  at  a 
number  of  RAS  trials-ranges  were  adapted  to  track  the  satellite*,  and  a  Pegasus 
computer  program  was  hurriedly  produced  for  the  analysis  [l]  of  observations  from  a 
single  transit.  Though  kinetheodolites  observe  azimuth  and  elevation,  the  program  was 
•written  [2]  to  work  with  direction  cosines  (to  which  conversion  is  trivial),  as  the 
original  objective  was  the  processing  of  data  from  a  radio  interferometer  specially 
built  at  Las ham. 

The  program  was  inevitably  somewhat  primitive  and  determined  only  four  orbital 
parameters,  quantities  equivalent  to  the  standard  elliptic  parameters  i,  Q,  w  and 
M  but  defined  in  relation  to  a  crossing  (by  the  satellite)  of  the  latitude  of  the 
observing  station.  This  restriction  was  necessary  for  program  stability,  and  was 
convenient  in  practice  since  the  remaining  two  parameters,  equivalent  to  a  and  e  , 
could  be  obtained  from  visual  observation,  over  a  number  of  transits,  of  the  orbital 
period,  together  with  an  assumption  about  the  slow  variation  of  the  perigee  radius. 

Results  of  geophysical  significance  were  obtained  when  the  Pegasus  program  was 
used  to  analyse  kinetheodolite  observations  of  Sputnik  2,  results  that  were  first 
recorded  in  Refs  1  and  2.  First,  the  (secular)  variation  in  ft  ,  over  a  period  of 
months,  was  0.7%  less  than  had  been  predicted;  this  result,  which  was  controversial 
at  the  time,  was  obtained  by  fitting  a  value  of  Q  to  all  the  individual  values  of 
0  derived,  and  it  amounted  (in  today's  terminology)  to  a  first  equation  of  condition 
on  the  even  zonal  harmonics  of  the  Earth's  gravitational  field.  Secondly,  an 
unexpected  secular  variation  in  i  was  detected,  amounting  to  a  daily  decrease  of 
about  10-3  deg/d;  this  variation  was  later  interpreted  as  due  to  the  rotation  of  the 
atmosphere  and  its  detection  was  the  starting  point  for  the  subsequent  RAE  discovery 
that  the  rotation  is  (at  least  for  the  lower  'upper  atmosphere')  actually  a  super¬ 
rotation.  These  and  other  geophysical  results  from  the  early  satellites  have 
recently  been  reviewed  by  King-Hele  [3]« 


•  The  author  was  then  working  at  RAE  Orfordness,  where  the  majority  of  the  kinetheodo¬ 
lite  tracking  was  undertaken,  and  took  part  in  the  first  tracking  of  Sputnik  2  at 
about  5  am  on  9  November  1957  (six  days  after  launch).  He  was  responsible  for  super¬ 
vising  the  film  reading  and  forwarding  the  data  (angles  of  azimuth  and  elevation)  to 
the  parent  establishment  at  Farnborough.  Early  in  1958  he  transferred  to  Farnborough 
(GW  Department)  and  became  responsible  for  the  preliminary  processing  of  much  of  the 
data  from  Orfordness  and  the  other  RAE  outs tat ions,  for  use  in  Merson's  original 
computer  program. 


THE  PEGASUS  PROGRAM  BASED  ON  MERSON 'S  PERTURBATION  THEORY 


Though  striking  successes  were  achieved  by  the  program  that  has  just  been 
described,  perturbations  were  entirely  neglected.  This  was  permissible  because  each 
orbit  determination  used  observations  from  only  a  very  short  arc  of  the  orbit;  but 
by  the  same  token,  the  kinetheodolite  observations  were  not  really  being  exploited 
to  their  full  accuracy  -  the  accuracy  of  the  derived  orbital  parameters  was 
equivalent  to  a  positional  error  of  between  1  and  10  km.  Thus  the  inherent  limitations 
of  the  early  program  were  too  deep-rooted  for  more  than  interim  use  -  ie,  until  such 
time  as  a  proper  perturbation  theory  could  be  developed,  and  then  implemented  on  the 
Pegasus  computer. 

An  appropriate  theory  was  developed  by  Merson  [4]  in  i960  and  implemented  during 
1961  by  Merson,  Gooding  and  Tayler[5,6]. 

The  orbital  elements  in  this  theory  are  neither  osculating  nor  mean  (in  the 
commonly  understood  sense)  but  'smoothed',  the  significance  of  this  being  that  first- 
order  short-period  perturbations  due  to  J 2  can  be  fully  represented  by  a  variation  of 
smoothed  elements  that  is  greatly  reduced  in  comparison  with  the  variation  of  osculating 
elements  -  it  entirely  vanishes  for  circular  orbits.  Because  the  short-period  variation 
is  not  entirely  eliminated,  as  it  is  when  mean  elements  are  used,  it  was  decided  that 
smoothed  elements  should  only  be  published  at  ascending  nodes  of  the  orbit;  thus  'time 
at  a  node',  rather  than  'mean  anomaly',  is  the  sixth  orbital  element  (the  others  being 
a,  e,  i,  0  and  <u  )  in  the  computer  program,  and  elements  were  determined  at  regular 
nodal  intervals,  usually  every  25th  or  every  50th,  rather  than  at  epochs  uniformly 
separated  in  time.  The  theory  ignores  short-period  perturbations  due  to  ,  with 
l  >  2  ,  but  Ref  4  gives  explicit  formulae  for  the  secular  and  long-period  perturbations 
due  to  Jy  J^,  and  Jg  ;  corresponding  formulae  are  given  for  the  second-order 

effects  due  to  J.,  ,  both  in  osculating  and  smoothed  elements  -  the  effect  on 
smoothed  e  was  wrongly  given  as  zero  in  Ref  4  and  this  was  corrected  in  a  subsequent 
RAE  version  of  the  paper  (Technical  Note  Space  26).  In  implementing  the  theory, 
however,  no  distinction  was  made  between  secular  and  long-period  perturbations:  both 
types  of  perturbations  had  to  be  represented  via  the  coefficients  of  (optional)  poly¬ 
nomials  (of  degree  up  to  5)  that  were  provided  for  e,  i,  fi  ,  and  the  mean 
motion  n  ;  the  'constants'  in  these  polynomial  were  (except  for  n  )  the  values  of  the 
smoothed  elements  at  the  epochal  node,  and  the  coefficients  could  be  either  determined 
by  the  program's  differential -correction  component  or  held  fixed  at  predetermined 
values.  In  practice,  first-degree  coefficients  for  e,  i,  fl  and  u>  were  normally 
held  fixed  at  values  derived  from  the  prior  running  [6]  of  separate  programs  written 
for  the  purpose  (one  to  give  the  contributions  to  these  coefficients  due  to  and  J  , 
one  to  cover  lunisolar  effects,  and  one  to  cover  certain  effects  of  atmospheric  drag 


and  rotation);  occasionally,  second-degree  coefficients  would  be  used  in  the  same  way. 
Since  Merson's  theory  dealt  only  with  effects  due  to  the  Earth's  gravitational  field, 
lunisolar  effects  were  based  on  formulae  given  by  Cook  [7]  and  Gooding  [3].  Drag 
effects  were  represented  in  combination  between  the  main  program  and  the  separate 
program:  in  the  main  program,  the  along-track  acceleration  could  be  represented 
empirically  by  the  coefficients  of  the  polynomial  for  n  ,  whilst  the  variation  in  e 
would  be  constrained  (apart  from  any  polynomial  effect)  to  make  a(1  -  e)  constant;  the 
separate  program  gave  polynomial-coefficient  values  for  the  variation  of  e  and  to 
also  i  and  ft  to  cover  atmospheric  rotation  -  proportional  to  the  n  variation  on 
the  basis  of  an  assumed  value  for  density  scale  height. 

The  program  described  was  used  for  orbit  determination  for  a  number  of  satellites 
between  1962  and  1968,  a  good  example  being  for  Ariel  2  (196^-15A)  [9].  Its  main 
limitation  was  that  only  observations  expressible  as  right  ascension  and  declination 
could  be  accepted:  this  allowed  other  forms  of  directional  measurement,  in  particular 
azimuth  and  elevation  or  direction  cosines,  since  these  could  be  transformed  to  right 
ascension  and  declination,  but  the  inability  to  accept  range  data  was  a  serious  defect; 
there  was  no  possibility  of  range  rate  (Doppler)  measurements  being  used,  as  the  program 
did  not  even  compute  the  velocity  of  the  satellite,  the  theory  only  being  taken  as  far 
as  the  generation  of  position.  The  accuracy  was  much  superior  to  that  of  the  earlier 
program,  however:  typically,  orbit  determination  (and  hence  the  generation  of  inter¬ 
polated  satellite  position)  was  accurate  to  %  km  over  a  two-day  period  [9].  For  close- 
earth  satellites  the  limitation  on  accuracy  was  always  the  unpredictability  of  drag,  and 
it  was  this  -  not  the  failure  to  distinguish  between  secular  and  long-period  perturba¬ 
tions  -  that  restricted  the  applicability  to  periods  sometimes  as  short  as  two  days. 

The  inadequacy  of  the  long-period  representation  over  longer  periods,  as  well  as  the 
unsatisfactory  business  of  having  to  run  separate  programs  before  the  main  programs,  was 
soon  recognized,  however,  and  it  was  a  prime  candidate  for  improvement  when  (as  we  see 
in  the  next  section)  it  became  inevitable  that  a  new  program  would  have  to  be  written. 

Two  concluding  points  are  worth  making  about  Merson's  theory,  as  implemented  in 
the  Pegasus  computer  program.  The  first  is  that,  though  the  J 2  theory  was,  in  itself, 
free  of  singularity,  sources  of  singularity  were  present  in  the  program:  thus,  'time  at 
a  node'  is  an  unsatisfactory  element  for  an  equatorial  orbit,  where  nodes  are  undefined, 
and  the  same  is  true  of  the  elements  ft  and  «  ;  for  ft  and  o>  the  only  effect  is  an 
ill  conditioning  of  partial -derivative  matrices,  but  failure  in  nodal-time  definition 
is  more  serious,  as  this  undermines  epoch  itself.  Singularities  could  also  occur  in 
the  programs  run  before  the  main  program,  since  infinities  occur  in  the  rate  of  change 
of  ft  for  equatorial  orbits  and  in  the  rate  of  change  of  a>  for  circular  and/or 
equatorial  orbits. 

The  second  point  is  that  the  Pegasus  elements  were  in  certain  respects  superior 
to  the  elements  used  in  the  replacement  program.  As  remarked  in  Appendix  C  of  Ref  10, 


the  residual  short-period  variation  of  the  smoothed  elements,  which  was  a  nuisance,  could 

have  been  dealt  with  by  defining  a  set  of  'mean  smoothed  elements',  including  M  ,  to 

coincide  with  smoothed  elements  at  ascending  nodes.  They  would  then  qualify  as  'mean 

elements'  in  the  definitive  sense  of  a  recent  general  study  [ll]  and  it  is  possible  to 

state  the  values  of  the  k-constants,  in  equations  ( 'l 38)— < ^ ^ )  of  Ref  11,  that  are 

associated;  thus  k  ,  k  ,  k. ,  k_,  k  and  k„  are  given  by 

a  e  i  Q  cu  n 

k  =  ih  (2  +  3e2)  *  \  f  e2  cos  2uj  , 


ke  =  5h(5  +  2e2)  +  5f  008  2w  » 


k^  =  1  +  —  e  cos  «  , 


k_  =  2(v  -  M)  n  ~  e  sin  <u  , 
Q  u=0  ;>  1 


koj  =  -  M)u_q  +  j  J^32  eg  sin  w  -  (9f  +  4e2h)  sin  2  tuj 


k^  =  -3f  sin  2  <o  ; 


here  f  =  sin  i,  h  =  1  -  ^f  and  g  =  1  -  ^f  ,  whilst  the  suffix  in  (v  -  M)u_Q 
indicates  an  evaluation  at  the  ascending  node,  u  being  the  argument  of  latitude 
(a  »  ♦  v)  .  If  mean  smoothed  elements  had  been  defined  to  coincide  with  smoothed 
elements  at  descending  nodes,  rather  than  ascending  nodes,  the  signs  of  the  terms  in 
sin  <e  and  cos  w  would  be  reversed  and  the  suffix  on  (v  -  M)  would  be  'u=ir'  . 

It  then  follows  that,  if  the  elements  had  been  defined  by  'splitting  the  difference' 
between  ascending  and  descending  nodal  bases,  that  (in  particular)  k^  would  be  unity 
and  mean  inclination  would  have  coincided  with  the  most  natural  value  recommended  in 


■  1  -!f 


Ref  11}  ka  is  in  any  vase  the  recommended  value  (postscript  to  Ref  11),  so  two 
advantages  of  the  Pegasus  elements  (over  the  subsequent  PROP  elements)  have  been 
demonstrated.  A  full  comparison  of 'Pegasus  elements'  with  those  introduced  by  Kozai  [12] 
(and  hence  underlying  PROP)  was  made  in  Ref  13* 

3  THE  FORTRAN  PROGRAM  PROP 

There  were  three  overriding  reasons  for  wishing  to  replace  the  Pegasus  program. 
First,  being  written  almost  entirely  in  Pegasus  machine  orders,  it  could  not  survive  the 
eventual  demise  of  Pegasus  itself.  Second,  it  was  inconvenient  to  continue  using  (mean 
smoothed)  orbital  elements  that  were  not  recognized  outside  RAE.  The  SAO  (Smithsonian 
Astrophysical  Observatory)  were  using  a  computer  program  known  as  [l^]  001-3  (for 
differential  orbit  improvement),  based  on  Kozai ' s  theory  [12],  and  ESRO  (European  Space 
Research  Organisation)  were  in  the  process  of  adopting  [15]  the  same  theory  for  their 
DCQP  fl6l  (differential  correction  orbit  prograane).  When  a  new  RAE  program,  in 


Standard  Fortran  [17]  for  maximum  portability,  was  planned,  therefore,  in  1965,  it 
seemed  right  to  base  the  underlying  theory  (for  short-period  perturbations  due  to  J?  ) 
on  the  SAO/Kozai  mean  orbital  elements;  with  hindsight,  however,  this  was  a  retrograde 
step  in  relation  to  at  least  two  of  the  orbital  elements  ( a  and  i  ) ,  as  already 
mentioned  -  the  work  described  in  section  5  effectively  involves  a  return  to  Merson's 
a  and  i  .  The  third  overriding  reason  for  replacing  the  Pegasus  program  related  to 
its  main  limitation  in  use,  viz  the  restriction  to  angular  observations. 

The  proposed  program  was  given  the  name  PROP  (ora gram  for  the  refinement  of 
orbital  parameters)  and  its  dynamic  model  was  developed,  and  described  by  Merson  [18]. 

A  users'  manual  [lO]  was  in  due  course  written  for  the  PR0P3  version  of  the  program  and 
an  update  generated  for  the  PR0P6  version  [19 i20] .  The  program  has  been  extensively 
used  both  inside  and  outside  the  RAE  •«  it  has  been  estimated  that  some  10,000  individual 
orbit  determinations  have  been  carried  out  -  and  an  exhaustive  documentation  is 
available  [21 J  for  the  version  of  the  program  used  at  the  University  of  Aston  in 
Birmingham. 

An  important  feature  of  the  underlying  theory,  as  described  in  Hef  18 ,  is  that 
(like  the  Pegasus  program)  it  is  entirely  free  of  singularity,  even  though  the  normal 
elliptic  elements  (a,  e,  i,  fl,  to  and  M  )  are  used.  This  was  achieved  by  arranging 
that  both  short-period  and  long-period  perturbations  would  be  applied  in  composite 
manner,  the  six  perturbations  actually  computed  being  5i,  5Q  sin  i,  ou  +•  50  cos  i, 
5p/2p,  6r/p  and  6r,  p  being  the  parameter  (semi-latus  rectum)  given  by 
p  =  a(1  -  e^)  .  This  was  a  significant  improvement  on  Kozai's  approach,  since  the  only 
composition  of  perturbations  in  Ref  12  is  of  5a,  5e,  5w  and  5M  into  6r  and  5u, 
and  the  long-period  6u  is  singular  for  equatorial  orbits;  further,  the  retention  of 
six  perturbed  quantities  in  PROP  permits  velocity  to  be  covered  as  well  as  position. 

As  the  program  developed,  however,  it  was  found  [20]  desirable  to  replace  the  last  three 
composite  perturbations  (long-period  only)  by  5a,  6e  and  e  5v  .  (An  altogether  more 
natural  approach  to  the  composition  of  short-period  perturbations  has  been  developed  in 
the  last  few  years,  as  part  of  the  post-PROP  studies  described  in  section  5») 

As  with  the  Pegasus  program,  short-period  perturbations  are  present  in  PROP  just 
for  the  representation  of  first-order  effects,  and  this  immediately  limits  the 

accuracy,  for  close-earth  satellites,  to  around  50  m.  Secular  and  long-period  perturba- 
tions  are  present  for  the  representation  of  J 2  effects  and  (first-order)  Jj  effects 
for  l  ranging  from  3  to  as  high  as  it  is  desired  to  go  -  the  current  PR0P6  version  of 
the  program  would  have  to  be  recompiled  if  it  were  desired  that  l  exceed  40.  The 
first-order  (secular  and  long-period)  effects  of  are  represented  completely  -  see 

Ref  11  for  a  discussion  of  'order'  in  this  context  -  and  this  means,  in  particular,  that 
the  representation  is  correct  for  orbits  that  are  near-circular  [22]  or  near-equatorial 
or  both.  The  representation  is  via  recurrence  relations  [8*18],  and  these  are  very 
efficiently  implemented  -  the  appropriate  quantities  (functions  of  the  epoch  elements 
a,  e,  i  and  at  )  are  built  up  at  the  beginning  of  each  iteration  of  the  differential- 
correction  process  and  there  is  no  unnecessary  recomputing. 


The  correct  representation  of  long-period  perturbations  is  an  important  feature  of 
:0P,  but  the  most  important  new  feature,  in  relation  to  the  preceding  Pegasus  program, 
is  been,  as  already  remarked,  the  ability  to  handle  observations  other  than  -angular: 
ius  the  program's  first  application  was  to  range  and  range-rate  observations  of 
;o  2  (1965-31A)  [23].  The  facility  to  process  no  less  than  16  types  of  observation  was 
tilt  into  the  original  version  of  the  program;  for  example,  Types  1,  2,  3,  ^  and  7  refer, 
•spectively,  to  observations  of  p  (range)  alone,  a  and  5  (right  ascension  and 
•clination),  1  and  m  (direction  cosines,  of  the  sort  covered  by  the  previous  programs), 
alone,  and  simultaneous  observations  of  p ,  a  and  5  .  Observations  of  azimuth  and 
Levation  are  not  processed  separately  but,  as  with  the  Pegasus  program,  converted  to 
and  6  at  the  input  stage,  corrections  for  refraction  etc  being  made  at  the  same 
Lme,  as  appropriate.  Input,  which  is  part  of  PROP  and  not  ancillary  as  with  the 
egasus  program,  is  associated  with  the  sources  of  the  observations  (ie  the  organizations 
oncerned)  rather  than  with  the  observation  types,  so  that  their  observation  formats  can 
e  used  directly.  Some  of  the  effects  of  this  dissociation  are :  that  there  are  several 
ources  of,  for  example,  Type  2  observations ,  all  distinction  having  been  lost  by  the 
tart  of  the  differential-correction  process;  that  for  some  of  the  types  of  observation, 
g  Tyne  5  (  &  and  5  ),  there  have  so  far  been  no  sources  of  observations,  so  that  these 
ypes  have  been  redundant;  but  that  new  sources  of  data  can  easily  be  accommodated  ('with 
ctivation  of  a  redundant  type  if  necessary). 

The  basic  orbital  elements  of  PROP  are,  as  already  intimated,  the  Kozai  meah 
laments  e,  i,  Q,  co  and  M  ,  the  element  a  being  defined  from  n  (=  M)  by  a 
2~modified  version  of  Kepler's  third  law.  It  is  important  to  realize  that  e 
similarly  i  ,  etc)  differs  from  osculating  e  only  in  respect  of  the  short-period 
erturbation;  the  use  of  a  'mean  mean'  element,  obtained  by  removal  of  a  nominal  long- 
eriod  perturbation  from  the  mean  element,  was  eschewed  from  the  start  so  that  the 
rogram  would  be  free  from  any  singularity  in  the  vicinity  of  the  critical  inclinations, 
nother  way  of  saying  this  is  that  long-period  perturbations  in  PROP  are  represented  by 
efinite  integrals,  the  lower  limit  of  the  integration  being  the  epoch  at  which  the 
lements  are  being  determined.  Epochs  in  PROP  are  constrained  to  be  midnights,  but  it 
ould  be  easy  to  remove  this  restriction. 

The  mean  elements,  e  etc,  are  allowed  to  vary  in  two  ways  as  t  departs  from 
ts  epoch  value,  tQ  .  First,  explicit  polynomial  variation  (of  degree  up  to  5,  with 
he  variation  of  M  inducing  a  derived  variation  in  n  )  is  allowed,  the  polynomial 
oefficients  being  denoted  by  e^  etc  (  j  =  0,  ...  ,  5)  •  Thus  eQ,  iQ,  Qq,  cu0  and 
are  the  values  of  the  mean  elements  at  epoch;  is  an  obligatory  sixth  element 

equal  to  aQ  and  thereby  defining  aQ  );  M^,  ...  ,  are  optional  elements  that 


in  practice)  permit  an  empirical  representation  of  the  along-track  acceleration  due  to 

tmospheric  drag;  components  of  e.  and  i,  are  determined  from  n.  (=  (j  +  1)  M.  ) 

j  j  j  J 

o  represent  residual  atmospheric  effects;  and  components  of  .  id  <0^  are  deter- 

ined  from  the  secular  effects  of  the  even  (including  the  effect).  Residual 

opponents  of  e^,  i^,  QJ(  and  a»^  ,  known  [10]  as  'exclusive  elements',  are  permitted, 

or  use  in  modelling  perturbations  that  are  not  explicitly  covered  by  PROP,  but  they  are 


almost  never  used  in  practice.  The  other  way  in  which  the  mean  elements  are  allowed  to 
vary  is  by  incorporation  of  the  long-period  perturbations,  these  being  defined  (as 
already  made  clear)  so  as  to  be  zero  at  epoch;  this  variation  is  implicit,  rather  than 
explicit,  however,  as  the  perturbations  are  applied  in  the  non-singular  maimer  that  has 
been  described.  (A  source  of  confusion  in  notation  arises  here:  mean  elements  are 
often  [ll]  denoted  by  e"  etc,  PROP'S  eQ  being  the  Kozai  e  at  epoch;  since  the  long- 
period  variation  in  PROP  elements  is  not  explicitly  added  to  the  polynomial  variation, 
however,  Refs  18  and  10  reserved  the  notation  e  for  the  polynomial  variation  -  the 
notation  has  been  avoided  here,  to  prevent  misunderstanding.) 

A  complete  summary  of  the  orbited,  model  of  PROP  is  given  in  Appendix  B  of  Ref  24, 
the  users'  manual  for  PREP  (program  for  realizing  enhemeris  nrint-outs),  which  is  the 
program  for  ephemeris  generation  paralleling  PROP  for  orbit  determination.  Other 
programs  based  on  the  PROP  model  are  referenced  in  Ref  20. 

For  each  basic  element  the  set  of  k  selected  elements  in  the  model  (polynomial 

coefficients  e.  ,  with  j  =  0,  ...  ,  k-1  )  are  divided  into  two  subsets:  the  first  m 
3 

coefficients  (0  S  m  $  k  ,  so  either  subset  can  be  empty)  are  parameters  of  the 

differential-correction  process,  whilst  the  remaining  (k  -  m)  are  held  fixed.  Partial 

derivatives  of  all  the  observations  are  computed  with  respect  to  the  parameters  (elements 

of  the  first  subset)  during  each  iteration  of  the  process.  The  computation  is  purely 

analytical,  allowance  being  made  for  first-order  secular  variation  due  to  J 2  ,  but  for 

no  other  source  of  perturbation.  At  two  critical  points  of  the  process  (during  each 

iteration),  PROP  can,  at  the  user's  option,  modify  the  standard  set  of  elliptic 

elements  so  as  to  avoid  convergence  problems  associated  with  near-circular  or  near- 

equatorial  orbits.  Full  details  are  given  in  Ref  10;  for  an  orbit  with  e  close  to 

zero,  one  modification,  which  avoids  an  ill-conditioned  derivative  matrix,  effectively 

involves  the  replacement  of  partial  derivatives  wrt  uiq  and  Mq  by  derivatives  wrt  u>  q 

ar.d  U  (=  M  +  «  ),  whilst  another,  which  avoids  the  necessity  for  non-linear  oara- 

meter  correction,  involves  the  temporary  introduction  of  elements  eQ  cos  o)q  and 

e  sin  «  to  replace  e  and  w  .  The  former  modification  relates  to  the  high 
0  0  0  o 

negative  correlation  between  <u  and  M  when  e  is  small;  the  PROP  user  should  also 

0  0  o 

respond  to  this  situation  when  tabulating  orbital  elements  obtained  by  the  program, 
since  the  prooer  procedure  is  to  quote  only  to  the  number  of  figures  that  are 

actually  meaningful,  while  quoting  the  composite  parameter  (Mq  +  u)q)  to  more  decimal 
places  than  would  be  appropriate  for  Mq  . 

PROP  is  still  in  regular  use.  Its  main  limitation  is  probably  the  total  absence 
of  any  explicit  representation  of  lunisolar  perturbations,  though  some  implicit  allowance 
can  be  made,  as  in  the  Pegasus  program,  via  the  'exclusive'  polynomial  coefficients  that 
have  been  described.  For  close-earth  satellites,  over  a  period  as  short  as  one  day, 
the  error  is  no  more  serious  (say  up  to  50  m)  than  that  due  to  the  neglect  of  short- 
period  perturbations  for  Jj  with  \  >  2  ,  but  the  effect  rapidly  becomes  more  serious 
for  more  distant  orbits. 

The  only  tesseral  (or  rather  sectorial)  harmonic  represented  is  2  ’  ®ven 
the  Pegasus  program  it  was  realized,  during  the  Ariel  2  analysis  [9],  that  the  dominant 


g-track  effect  of  the  J_  ,  perturbation  was  important  enough  to  be  added  to  the 
tal  model,  and  this  effect  was  built  into  the  earliest  version  of  FRCP.  In  due  course 

as  decided  that  the  dominant  (semi-diurnal)  J,  _  effect  in  all  the  elements  should 

d  )d 

epresented,  and  Ref  20  gives  the  details. 

The  overall  accuracy  of  which  PROP  is  capable  is  of  the  order  of  50  m.  The  accuracy 
distribution  of  observations  are  not  always  good  enough  to  achieve  this,  of  course, 
for  Ariel  4  (1979-109A),  which  was  a  good  example  of  the  definitive  use  of  FRCP  for 
.t  determination  to  be  associated  with  the  experimenters'  telemetry  analysis,  an 
iracy  of  120  m  was  quoted  [25],  the  main  source  of  error  being  the  Minitrack  data. 

A  final  point  to  mention,  concerning  PROP  in  relation  to  the  preceding  Pegasus 
jram,  is  that  a  change  was  made  in  the  choice  of  equinox  to  define  the  x-axis. 

lgh  both  programs  took  the  equator  of  date  as  xy-plane,  the  apparently  natural  choice 

equinox  of  date,  made  for  the  Pegasus  program,  was  abandoned  for  PROP  in  favour  of  a 
si-equinox  defined  by  projection  of  the  mean  equinox  of  1950.0.  The  basis  for  this 
rid  system  is  given  in  Appendix  E  of  Ref  20;  it  follows  Kozai  and  the  SAO. 

PROD.  A  COMPANION  PROGRAM  TO  PROP 

The  function  of  PROP  is  to  refine  a  set  of  Kozai  mean  orbital  elements,  and  usually 

e  other  parameters  to  represent  drag  effects,  by  fitting  to  observations  of  a  satellite 

interest,  collected  over  a  number  of  days.  Over  such  a  short  period  the  absence  of 
i3olar  perturbations  is  not  serious;  but  when  it  comes  to  interpreting  sets  of 
ments  (usually  at  uniform  intervals  of  time)  spanning  a  period  of  perhaps  several 
ths,  the  neglect  of  these  perturbations  can  in  most  cases  no  longer  be  tolerated, 
t  was  wanted,  therefore,  as  a  companion-  to  PROP,  was  a  program  that  v/ould  generate 
evolution  of  a  3et  of  mean  elements  over  long  periods,  allowing  for  lunisolar 
vitational  attraction  as  well  as  for  the  zonal  harmonics  of  the  Earth.  (Except  for 
its  experiencing  resonance,  referred  to  again  at  the  end  of  this  section,  or  when 
reme  accuracy  is  required,  it  is  legitimate  to  neglect  the  tesseral  harmonics,  since 
general  they  cause  perturbations  that  have  periods  of  less  than  a  day.)  PROD 
ogram  for  orbit  development)  was  developed  with  this  rationale. 

The  program  operates  by  numerical  integration  of  Lagrange's  planetary  equations, 

'  ariation  of  parameters'  method.  The  method  of  integration  is  a  fourth-order 
ge-Kutta  process,  which  is  not  very  efficient  over  long  periods  of  time.  No 
ficulty  has  been  found  in  practice,  however,  mainly  because  a  fairly  long  integration 
p  time  can  be  used  since  no  short -period  terms  are  carried  in  Lagrange's  equations. 

A  full  account  of  the  theory  underlying  PROD  was  given  by  Cook  [26],  and  the 
gram  has  been  regularly  used  for  assessing  the  lifetimes  (via  a  determination  of  the 
Igee-height  variation)  of  newly-launched  satellites  in  highly-eccentric  orbits.  The 
a  features  of  the  theory  and  program  are  now  outlined  together. 

T3ia  program  only  covers  the  elements  a,  e,  i,  Q  and  u  .  The  sixth  element  (M) 
i  included  in  the  theory,  and  it  was  intended  to  add  H  to  the  program,  but  in  practice 


:he  evolution  of  M  has  never  been  required  -  since  tesseral-harmonic  perturbations  are 
lot  represented,  the  rates  of  change  of  the  mean  elements  a,  e,  i,  1  and  '-u  are 
independent  of  M  .  (Further,  the  variation  in  a  is  always  set  to  zero.) 

Only  gravitational  perturbations  are  represented,  though  the  original  intention  was 

to  cover  drag  and  solar  radiation  pressure  as  well.  The  Earth's  zonal  harmonics  are 

covered  essentially  as  in  PROP,  with  first-order  secular  and  long-period  perturbations 

represented  uo  to  J .  ,  as  reouired,  together  with  second-order  J..  rerturbacions 
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there  is  an  error  in  the  rate  of  change  of  i  ,  as  quoted  in  Re;  26,  since  this 

rate  of  change  is  proportional  to  sin  2(u  not  cos  2<j  .  The  contributions  of  the  J 
to  the  rates  of  change  are  generated  by  recurrence  relations,  similar  to  those  in  PROF 
though  computed  somewhat  less  efficiently.  Since  the  standard  elliptic  elements  are 
integrated  directly,  there  is  no  avoidance  of  singularity  as  there  is  in  PROP  and  this 
is  probably  the  main  defect  of  the  program.  (A  distinct  program,  known  as  PR0D2,  was 
developed  to  cope  with  near-circular  and  near-equatorial  orbits,  but  this  is  too  slow 
for  general  use  as  the  Fortran  coding  makes  liberal  use  of  complex  variables.) 

Lunisolar  perturbations  are  represented  in  such  a  way  that  the  user,  via  four 
digits  of  one  of  the  quantities  on  the  control  card  (part  of  the  data  input),  has  a  pair 
of  options  for  both  the  lunar  and  the  solar  components  of  the  element  rates  of  change. 

The  first  option  relates  to  the  number  of  harmonic  terms  used  in  the  standard  expansion 
of  the  disturbing  function  due  to  the  attracting  body  -  see  equation  (10)  of  Ref  26  or 
equation  (22)  of  Ref  8.  A  single  term  should  always  be  adequate  for  the  solar 
attraction,  in  view  of  the  Sun's  great  distance,  but  at  least  two  will  normally  be 
required  for  the  lunar  attraction  (the  second,  involving  a  harmonic  function  of  degree  3, 
being  known  as  the  parallactic  term)  and  sometimes  (for  satellites  in  highly  eccentric 
orbits)  it  may  be  necessary  to  take  as  many  as  four  -  the  limit  allowed  is  eight,  ie 
harmond.es  up  to  degree  9.  The  other  option  permits  a  simple  choice  between  a  singly 
averaged  disturbing  function,  ie  averaged  only  with  respect  to  the  mean  anomaly  of  the 
satellite  (to  avoid  short-period  perturbations),  and  a  doubly  averaged  disturbing 
function,  ie  averaged  also  with  respect  to  the  mean  anomaly  of  the  attracting  body 
itself.  It  is  important  to  realize  that  it  is  not  just  a  matter  of  excluding  certain 
terms  when  double  averaging  is  opted  for,  since  there  is  a  fundamental  difference  in  the 
way  the  terms  of  the  disturbing  function  are  developed  -  with  single  averaging  the 
attracting  body  is  effectively  at  an  arbitrary  fixed  point  in  space  (to  be  specified  by 
a  subroutine)  at  each  integration  point,  whereas  with  double  averaging  it  is  only  the 
orbit  of  the  body  that  is  relevant,  the  body  being  effectively  replaced  by  a  ring  of 
mass  distributed  around  that  orbit.  Far  from  requiring  fewer  terms,  this  means  that 
double  averaging  actually  requires  more;  thus  equations  (59)-(63)  of  Ref  20  involve 
triple  summation,  whereas  equations  (3"0-(35)  involve  only  double  summation.  In 
practice,  because  the  Earth-Sun  orbital  period  is  so  long  (as  well  as  because  fewer 
terms  are  involvedl),  it  will  be  normal  to  take  the  solar  function  with  only  single 
averaging;  for  the  lunar  function,  on  the  other  hand,  it  will  often  make  sense  to  opt 
for  double  averaging,  since  it  would  otherwise  be  necessary  to  restrict  the  integration 
step  time  to  being  a  very  small  fraction  of  a  month. 


The  most  serious  limitation  of  PROD  is  the  absence  of  all  tesseral-harmonic 
representation,  which,  together  with  the  omission  of  sixth-element  variation,  precludes 
the  coverage  of  resonance  effects.  Such  effects  can,  of  course,  be  of  critical  import¬ 
ance,  in  particular  for  24h  and  12h  orbits.  It  is  worth  mentioning,  therefore,  another 
program  that  is  unrelated  to  either  PROP  or  PROD.  Known  as  TCSKEF,  it  originated  as  a 
program  used  by  the  Royal  Air  Force  for  the  quick  generation  of  ephemerides  for  Skynet 
communications  satellites,  but  was  taken  over  by  RAE  and  made  available  as  a  general- 
purpose  tool  for  predicting  mean  orbital  elements,  ephemerides  and  look  angles.  The 
orbit  generator  has  been  documented  elsewhere  £27]  and  it  suffices  here:  first,  to 
remark  that,  though  not  totally  analytic  like  PROP,  TCSKEF  is  much  more  analytic  chan 
PROD,  in  that  it  is  finite  increments  to  elements,  rather  than  their  rates  of  change, 
for  which  formulae  are  available;  and  secondly  to  summarize  the  sources  of  perturbation 
covered.  Only  one  zonal  harmonic,  ,  is  included,  and  this  of  course  is  an  immediate 
limitation  on  applicability;  as  described  in  Ref  27,  the  main  24h  resonance  effects  are 
covered  via  the  tesseral  harmonics  J_  _  ,  J,  ,  and  J,  .  ,  and  the  main  12h  effects 

via  J  .  ,  but  the  orbit  generator  has  since  been  enhanced  so  that  some  additional 
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effects  of  resonance  are  explicitly  covered  auid  others  can  easily  be  added;  lunisolar 
attraction  is  represented,  but  only  for  the  leading  harmonic  term  in  the  expansion  of 
the  disturbing  function,  by  a  procedure  that  gives  an  appropriate  weighting  to  the 
effects  associated  with  single  and  double  averaging;  a  primitive  representation  of  the 
effects  of  eclipse-free  solar  radiation  is  included. 

5  RECENT  DEVELOPMENTS 

After  PROP  and  PROD,  work  at  RAE  concentrated  on  orbit  generation  by  Cowell's  method, 
also  known  rather  misleadingly  as  'special  perturbations',  which  is  at  the  opposite  end  of 
the  spectrum  from  the  'general  perturbations'  method  underlying  PROP,  since  it  is  based  on 
pure  (and  highly  efficient)  numerical  integration  -  rates  of  change  of  position  and 
velocity  coordinates  are  integrated,  so  the  method  does  not  recognize  such  a  concept  as 
'orbital  element'.  Cowell's  method  is  employed  for  its  superior  accuracy,  and  the  pair  of 
programs  used  with  the  Skynet  2  communications  satellites  were,  for  orbit  determination  and 
ephemeris  generation  respectively,  SP0D2  [28]  and  SKEPH2  [29] •  Generalized  versions  of  these 
programs,  permitting,  eg,  the  use  of  arbitrarily  extensive  geopotential  models,  have  been 
developed  since  then;  known  as  POLO  and  PINTO  respectively,  they  have  not  yet  been  docu¬ 
mented.  No  more  will  be  said  about  any  of  the  programs,  however,  as  their  purely  numerical 
nature  puts  them  outside  the  scope  of  the  present  paper. 

It  occurred  to  the  author  early  in  1980,  after  completion  of  the  orbit  theory  that 
was  in  due  course  published  as  Ref  11,  that  a  new  approach  to  orbit  generation  (and 
eventually  orbit  determination)  was  possible,  which  would  combine  the  speed  and 
sophistication  of  analytical  generation  with  the  accuracy  of  pure  numerical  integration. 

A  mixed  analytical-numerical  procedure  would  thereby  be  adopted  that  could  be  described 
as  semi -analytical ,  though  (just  as  with  the  program  TCSKEF,  referred  to  at  the  end  of 
the  last  section),  the  numerical  element  of  the  procedure  would  be  relatively  slight, 
amounting  rarely  to  the  splitting  of  any  excessively  long  period  of  time,  over  which 
mean  elements  are  to  be  propagated,  into  subperiods. 


There  is  nothing  intrinsically  original  in  the  idea  of  a  mixed  procedure.  The 
novelty  of  the  approach  lies  in  the  particular  interface  c.  n  between  the  analytical 
and  numerical  components  of  the  orbit  generator.  Two  aspects  axe  of  fundamental  import¬ 
ance  here.  The  first  concerns  the  treatment  of  short -period  perturbations,  which  is 
purely  analytic  and  involves  their  expression  in  a  (geocentric)  system  of  spherical- 
polar  coordinates  based  on  an  optimally  defined  mean  orbital  plane.  (Ref  11  was  based 
on  cylindrical  coordinates,  but  when  the  complete  set  of  formulae  for  perturbations 

were  developed  [po]  it  was  realined  that  spherical  coordinates  were  better.)  The  second 
aspect  of  the  interface  concerns  mean  elements,  which  are  conceptually  defined  so  as  to 
be  free  of  terms  involving  the  satellite's  mean  anomaly  (the  normally  defined  short- 
period  perturbations)  and  of  terms  involving  the  Greenwich  sidereal  angle  (ir.  perturba¬ 
tions  due  to  the  tesseral  harmonics),  except  inasmuch  as  these  quantities  appear  together 
in  terms  that  must  be  regarded  as  resonant  -  an  arbitrary  criterion  for  resonance  is 
required  here.  The  precise  definition  of  the  mean  elements  involves  the  assignment  of 
arbitrary  constants  in  the  analytical  integration  (the  k-constants  of  section  2  and 
Ref  11,  as  far  as  perturbations  are  concerned),  and  there  are  significant  rewards 

(in  the  simplicity  of  the  expressions  for  short-period  perturbations)  if  the  right 
assignment  is  made  -  the  assignment  for  mean  i  and  Q  effectively  defines  the  mean 
orbital  plane. 

Results  from  a  preliminary  test  of  the  new  approach  were  reported  [jl]  at  the 
19S0  COSPAR  meeting  in  Budapest.  The  test  program  covered  the  potential  field 
associated  with  and  (as  specimen  zonal  harmonics),  2  ^as  a  Prototype 

tesseral  harmonic),  the  Sun  and  the  Moon.  Short-period  perturbations  (in  cylindrical 
coordinates)  were  restricted  to  e-independent  terms  and  the  program  was  not  singularity- 
free  for  equatorial  orbits.  The  results  were  promising,  however,  the  main  source  of 
error,  for  an  orbit  of  period  4b,  eccentricity  0.002  and  inclination  70°,  being  an 
unmodelled  coupling  between  o'2  and  J2  2  that  led  to  an  along-track  error  signature 
of  period  6  h  and  amplitude  about  0.7  m;  the  errors  were  much  greater  for  a  higher 
orbit  (l8h  period  for  the  test  case),  the  inadequately  modelled  factor  then  being  the 
motion  of  the  Moon. 

The  next  stage  of  the  work  was  to  complete  the  analysis  of  short-period  perturba¬ 
tions  of  second  order  in  J 2  and  first  order  in  ,  to  make  them  valid  for  (elliptic) 
orbits  of  arbitrary  eccentricity,  at  the  same  time  removing  the  source  of  singularity.  The 
results,  now  most  compactly  expressed  in  spherical  coordinates,  were  presented  [30]  at 
the  1982  IAF  meeting  in  Paris;  for  an  orbit  of  period  4  h,  eccentricity  0.5  and  inclina¬ 
tion  40°,  with  only  J2  and  represented,  the  dominant  error  (against  an  integrated 
yardstick  ephemeris)  was  found  to  be  in  apogee  position,  growing  at  about  6  cm  per 
revolution  and  due  to  coupling  between  J 2  and  .  A  full  account  of  the  underlying 
analysis  is  to  be  issued  shortly  [32]  and  work  is  now  concentrated  on  the  coverage  of 
J2  2  for  orbits  of  arbitrary  (<  1)  eccentricity. 

The  eventual  goal  is  a  fast  program  that  would  replace  the  PROP/PROD  combination 
and  provide  an  accuracy  of  order  1  ra.  Many  steps  remain  to  be  taken  before  this  goal 
can  be  seen  a a  realistic,  and  some  of  these  are  listed  aa  a  conclusion  to  this  paper. 
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(a)  The  first  order  coverage  of  J,  must  be  extended  to  an  arbitrary  zonal  harmonic; 
implementation  of  this  will  only  be  straightforward  if  recurrence  relations  can  be  found 
for  the  perturbations  in  spherical  coordinates. 

(a)  The  present  J.  .  study  must  be  brought  to  a  satisfactory  conclusion,  and  then 
extended  to  arbitrary  tesseral  harmonics,  again  with  the  use  of  recurrence  relations. 

(c)  A  sufficiently  accurate  representation  of  lunisolar  perturbations  must  be  included 
in  practice  it  is  likely  that  the  epithet  'sufficiently  accurate*  will  be  defined  on  the 
basis  of  an  upper  bound  on  a  and  e  for  the  specified  accuracy  to  be  achievable. 

(d)  Atmospheric  drag  must  be  covered,  unless  the  program  is  to  remain  unusable  for  the 
vast  number  of  orbits  below  some  (perigee)  altitude  that  is  itself  a  function  of  the 
satellite's  mass/shape  characteristics  and  the  solar-cycle  variation  in  atmospheric 
density. 

(e)  Other  perturbations,  such  as  those  associated  with  the  tides  and  those  due  to 
solar  radiation,  must  be  covered,  as  well  as  the  effects  of  the  Earth's  precession  and 
nutation. 

(f)  Sufficiently  accurate  partial  derivatives  must  be  programmed ,  so  that  the  program 
routinely  converges  to  the  right  results. 
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